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Abstract 

This paper is concerned with the motion of a time dependent hypersur- 
face d£l(t) in R d that evolves with a normal velocity 

V n = n — T n da, 
Jdn(t) 

where k is the mean curvature of dH,(t), and j T stands for -X, J Jm Phase field 
approximation of this motion leads to the nonlocal Allcn-Cahn equation 

d t u = Au- \W(u) + \ / W'{u) dx, 
e e Jq 

where Q is an open box of R d containing dil(t) for all t. We propose a 
modified version of this equation: 

d t u = Au-^W'(u) + ^ v / 2W(u)(^J y/2W(u) dx^j J W'(u)dx, 

and we show that it has better volume preserving properties than the clas- 
sical one, even in the presence of an additional forcing term g. 



1 Introduction and motivation 

In the last decades, a lot of work has been devoted to motions of interfaces, and 
particularly to motion by mean curvature. Applications concern image processing 
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(denoising, segmentation), material sciences (motion of grain boundaries in alloys, 
crystal growth), biology (modelling of vesicles and blood cells). 

In this paper, we are interested in phase field equations as an approximation 
to motion by mean curvature with a forcing term and a volume constraint. 

For t in [0, T], let O(t) denote the evolution by mean curvature with a forcing 
term of a smooth bounded domain Qo i n More precisely, the normal velocity 
V n , with normal n pointing towards the exterior of f2(t), is given at a point x of 
dn{t) by 

V n = K + g, (1) 

where k denotes the mean curvature at x, with the convention that k is negative 
if the set is convex, and where g = g(x, t) is a given smooth forcing term. In this 
work, we only consider smooth motions, which are well-defined if T is sufficiently 
small [2J. Singularities may develop in finite time, however, and one may need to 
consider evolutions in the sense of viscosity solutions [3] |8] . 

The evolution of O(t) is closely related to the minimization of the following 
energy: 

J(Q) = / Ida- [ gdx. 
Jdn Jn 

Indeed, one can view ([1]) as a first order optimality condition for this energy. The 
functional J can be approximated by a Ginzburg-Landau energy [IP], 19]: 

J e (u) = f ( -\\7u\ 2 -\ — W(u) J dx — cw f gudx, 
7rA 2 e J JR d 



where e is a small parameter, W a double well potential with wells at and 1, 
I, 

2' 



for example W(s) = \s 2 {l — s) 2 , and where 



cw = f yj2W{s) ds. 
Jo 

Modica and Mortola [101 [9] have shown the T-convergence of J e to cy/J in L 1 (R d ) 
in the absence of forcing terms (see also [3]). The extension of these results 
to motions with bounded forcing terms is straightforward. The corresponding 
Allen-Cahn equation obtained as the gradient flow of J e , reads 

d t u = Au- -^W'{u) + ^c w g. (2) 

This equation is usually solved in a fixed box Q of R d , which contains the motion 
Q(t) for all t in [0,T]. Existence, uniqueness and a comparison principle have 
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been established for this equation (see for example chapters 14 and 15 in [2]). To 
this equation, one usually associates the profile 



q = argmin j^Q7 /2 + W{-y)^j ds ; 7 6 v\ , 

where V is the space of functions in iJ^ c (R) that satisfies 7(— 00) = 1, 7(+oo) 
0, 7(0) = \. For t in [0,T], the motion can be approximated by that of 



a(t) 



|xGR d ; u £ (x,t) > ^J, 



where u e solves ([2]) with the initial condition 

' d(x, r^o) 



« e (x,0) 

Here d(x,Cl) denotes the signed distance of a point x to the set Q. The conver- 
gence of <90 e (i) to dd(t) has been proved for smooth motions [5] and in the 
general case without fattening [3l [8] . The rate of convergence has been proven to 
be 0(e 2 |log e| 2 ). Actually, a formal asymptotic expansion shows that u e behaves 
like 

„ e(l , t) ^(^hM))) +£9( ^*fM5)) +0(e »,, (3) 

where 77 is defined as the solution in H? CR), with polynomial growth, of 

U"-W"(q)ri = -c w + q', 

\v(o) = 0. 

When 5 = 1, the modified profile s 1— ► = g(s) + er/(s) (see figure [1]) can be 
evaluated at s = ±00, where it takes the respective values 

e TT ^„. . and 1 + 



'W(l) W"(0) 

These values correspond to the positions of the wells of a modified double well 
potential W £)ff , defined by W' eg = W — ecy/9 and W £iff (0) = 0. 

Our main interest is the numerical simulation of interfaces dQ(t) evolving 
from <9f2o with normal velocity given by 

V n = K + g- f {K + g)da. (5) 
JdClit) 
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Figure 1: Profile of the function s i— » q(s) + er](s) when W(s) = \s 2 (l — s) 2 . 

In this case, it is easy to see that the volume of f2(i), 

|n(t)| = f ltte, 
7n(t) 

remains constant in time. For instance, using the results in |12j . one may check 
that the shape derivative of the volume is zero. The usual strategy to approximate 
(J3J) is based on the remark that the mass 

/ u e dx 

jR, d 

is a good approximation of the volume One can then add to the Allen- 

Cahn equation an extra forcing term A(t), independent of x, in order to impose 
the conservation of mass. This leads to the following equation: 

dtu = Au ^ (W(u) — ecwg) H — cw^- 



The forcing term A can be viewed as a Lagrange multiplier associated to the 
volume constraint. It can be determined by integrating the equation over Q, 
which gives 



f ~2 ( w '( u ) ~ tcwg) dx. 



CW JQ 

In the case where g = 0, the previous equation reduces to 



d t u = Au - \w'(u) + f \w\u) dx, (6) 
6 JQ e 



which is the classical Allen-Cahn conserved equation (see [TT] and [6]). Formally, 
one can think of this equation as an approximation to motion by mean curvature 



with a modified forcing term g e (t), independent of x and given by 



e 



<? e = — -f — W (u) dx. 
c w JQ 

In view of expansion ([2D, one expects solutions of © to behave like 



By integration over Q, one sees that (see proposition Q] further) 
1 u e dx = \n e (t)\+eg e J V ^ X,n e e{t) ^ dx + 0(. 
As the mass of u e is conserved, as g t = 0(1), and as 



= 0(1), 

given the values of 77 at ±00, one expects that 

|n e (t)| = |n | + o(e) 

only. This is not satisfactory for many applications, where loss of volume during 
numerical computations strongly affects the dynamics. 

The aim of this work is to propose another phase field model that has better 
volume conservation properties than the conserved Allen-Cahn equation. The 
paper is organised as follow: 

In section [21 we introduce the following phase field approximation for mean 
curvature flow with a forcing term: 

d t u = ku-^(w\u)-e^/W(u)g). (7) 
It can be seen as the gradient flow of 

J e (u)= I ( -\Vu\ 2 + -W(u) J dx- [ G(u)gdx, 
JrA 2 e / JR, d 



with 



G(s) = [ y/2W(t) dt. 
Jo 



We first explain via a formal asymptotic analysis why solutions of ([7]) are expected 
to take the form 

Ul(x , t)= J d J^M>l) + 0^). (8) 
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Then, following an argument due to [5], we rigorously prove the convergence of 
this phase field equation to the motion (pQ). 

In section [3l we consider the evolution Q(t) of a smooth bounded domain f^o 
according to 

V n = K + g - f (K + g)da. 
J dn 

Let u f be a solution of 



d t u = Au - \ (W'(u) - ey / 2W(u)g 



e 2 



+ ~ 2 r / (w\u)-e^WW)9)d X . (9) 



We show that if u e behaves like in expansion ([8]), then for all t in [0, T], 

\tt \ = / u e (x,0)dx + O(e 2 ) 

= / u e (x,t)dx + 0(e 2 ) 
Jn d 

= \n e (t)\ + o(e 2 ), 



while the solutions of ([2]) may only conserve volume up to order e. 

In section^ we present numerical evidence for the above claims, which show 
that the modified phase field model ([9]) has indeed better volume preservation 
properties. 

2 A modified reaction— diffusion equation for mean 
curvature flow with a forcing term 

Let dQ(t) denote an evolving hypersurface of codimension 1 in R d , with velocity 
law V n = k + g. This motion can be interpreted as the energy gradient of 



J(U) = [ Ids - [ gdx. 
Jan Jn 



Let W be a bounded double well potential. In this whole section, we will 
for convenience use a potential with wells at —1 and 1, for example W(s) = 
min{^(l — s 2 ) 2 ,M}, where M is a given positive constant. Our strategy is to 
introduce a modified Ginzburg-Landau energy J e defined on L 1 (R d ) by 



f ( tt|Vu| 2 + -W(u) ) dx - [ G(u)gdx ifuetf 1 ^ 
+oo otherwise, 
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with 

G(s) = [ y / 2W{t) dt. 
Jo 

The function G is Lipschizt continuous. For g in L°°(R d ), the term u i— > 
J Rd G(u)g dx acts as a continuous perturbation in the L 1 (R ci ) topology of the 
classical Modica-Mortola energy. The stability of T-convergence with respect 
to continuous perturbations allows us to extend the Modica-Mortola result to 
the case at hand, and show that J t T-converges to cy/J- The gradient flow of J e 
should then provide a mean to approximate the motion of dfl(t) via the resolution 
of the reaction-diffusion equation ([7]) . 

Remark 1. In the simplest case where g = 1, equations ([2]) and (|7|) can be 

expressed as Allen-Cahn equations with particular double well potentials respec- 
tively equal to Wi >e (s) = W{s) + ecws, and W2,e(s) = W(s) + eG(s). These 
two potentials are related through the position and height of their wells, which 
are asymptotically equal as e — > 0. This explains why we expect that ([2]) and ([7]) 
converge to the same motion. 

2.1 Formal asymptotics for the modified Allen— Cahn equation 

We denote by u e the solution of equation ([7]) : 

8 t u = Au- ^W'(u) + -^2W{u)g, 

with initial condition 

u(x,0) = ~( d{xM 



e 

Our aim is to propose an asymptotic analysis of u e in the simplest two-dimensional 
radial case. Using polar coordinates (r,8), we consider a forcing term g which 
does not depend on 9: g = g(r,t). The initial set Qq is taken as a disk of radius 
1: 

n = {(r, 9) G [0, +oo) x [0, 2vr) ; r < 1} . 

Let Q,{t) be the mean curvature flow evolving from Oo according to the law 
V n = k + g. It is well known that in this case, £l(t) remains a circle for all t 
(recall that the forcing term g is supposed to be radial). We will denote by R(t) 
the radius of £l(t), solution of the following ODE: 

R' + ^=g(R,t), 

with initial condition i?(0) = 1. In this simple case, the solution u e is also radial 
and depends only on r. It satisfies 

d t u e - -d r (rd r u e ) + \w'(u t ) - -y/2W(u e )g = 0. 
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As u t is radial, every of its level sets is circular, and we denote by R t (t) the radius 
of {u e (r,t) = I}. We thus have R e (0) = 1 and u e (R e ,t) = |. We introduce the 
classical stretched variable y = r ~ £ (see [5]), and we define U e by 

U e (y,t) =u e (R e + ey,t). 

This new function U e satisfies 

d t U e - - e R' e d y U t - ^d y U e - ^d yy U e + (U e ) - 2W {U t ) g = 0. (10) 

We now consider asymptotic developments of U e and R t as follow: 

+00 +00 

us, t) = J2 *)> Re(t) = 

i=0 i=0 

with U o (0,t) = 5, Rq(0) = 12(0), and ?7<(0,t) = 0, 1^(0) = for all i > 1. We 
have 

W'(U e ) = W'(U ) + eW'iUop! + e 2 (V"(£/o)£/i + ^"(E/b)^) + 0(e 3 ), 

yW({7 ) 

g(r, t)=gley + Y, e*Ri,t j = g{R ,t) + ed r g(R ,t)(y + 12i) + 0(e 2 ). 



Using these equalities, (jlOp rewrites 
= l(5 ro C/ -W(C/o)) 



1 /„ „ „,„,„ 1 



- d t U + R[d y U + R' d y Ui + 5 ra C/ 2 + -^dy^i - - ^ ^ d y U ^-q 

+ g(R ,t) T^= U X + d r g(Ro,t)(y + R 1 )^2W{U ) 

- W'"(Uo)Ux - W"(U )U 2 
+ 0(e). 

Following powers of e, we will now identify each term to zero. 
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Terms in e . The first term Uq satisfies d yy Uo = W'{Uq) with initial condition 
U (0,t) = \ for all t in [0,T]. It can thus be identified to the profile q: 

VyGR, Vte[0,T], U {y, t) = q(y). 

Terms in e _1 . Knowing by definition of the profile that q' = —y / 2W(q), it 
follows from Uo(.,t) = q(.) that d y Uo = —\J2W{Uq). Equation ([IT]) then gives 

SyyUi - W'iUop! = -d y U [R' + -L - g(R ,t)^ . 

Multiplying this equality by d y Uo and integrating over R, we get 

R '° + 1T ~ 9{Ro,t) ) f^ d y u rf d y = ~ f ( d yv Ul ~ w"{u )u^d y u Q dy 



l dyfdyyUo-W'iUo^Uxdy 



= 0. 

As fn(d y Uo) 2 dy is strictly positive, we get the following equation on Rq: 

K + — =g{Ro,t), 

with initial condition Ro(0) = R(0). Hence -Ro can be identified to R since they 
both satisfy the same ODE with the same initial datum. It follows from (llip that 
U\ is solution of 

OyyUx ~ W"(U )Ut = 0. 

We then know (see section 3 in [5]) that there exists a(t) £ R such that Ui(., t) = 
a(t)q'(.). Indeed, the kernel of the operator A : iJ x (R) — » i2" _1 (R) defined by 
A( = C" — W"(q)C can be identified to span(g'). Using Ui(0,t) = 0, we conclude 
that a{t) = for all t, so that 

VyeR, Vte[0,T], Ui(y,t) = 0. 



Terms in e°. Using U\ = and dtUo = 0, we get from (llip that 

d yy U 2 - W"(U )U 2 = d y Uo (^TT^ " R i + 9 r9(Ro, t)(y + #i)) ■ 
Multiplying by d y Uo and integrating over R, we have 

R!i - - d r g{Ro,t)R^j J n (d y U ) 2 dy 
~~ jS dyyU2 ~ W "^ U2 ) d y Uod y + (^2 +9rg{Ro,t)^ J^y(d y Uo) 2 dy. 
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The first term in the right member vanishes as previously for U\. The second 
term also vanishes since d y Uo(.,t) = q'(.) is even with our choice of W. We 
deduce that R\ is solution of 



R 'i = (r^2+d r g(Ro,t)\Ri, 



with initial condition i?i(0) = 0. Hence R\(t) = for all t in [0, T]. Finally, U 2 
is obtained as the solution of 

d yy U 2 - W"(U )U 2 = yd y Uo ( + d r g(R ,t) 

\Ro 

Introducing the solution £ in ff 1 (R) of 

e(y)-w"{ q (y))ti(y) = yq'(y) 

with £(0) = 0, we can express U 2 : 

VyER, ViG[0,T], U 2 (y, t) = (-^ + d r g(Ro, t) 

We finally conclude from this formal asymptotic analysis that u e , solution of 
((7j), is expected of the form 

„ <(M) = ,(4^) +£2 (^ +8r9W , t) ) ? (^(£^) +0 (^ (12) 

where f2 e (£) converge to Q(t) in 0(e 2 ). 

2.2 Proof of convergence for the modified phase field model 

In this section, we closely follow the work of [Sj to prove the following theorem: 

Theorem 1. Let £l(t) be a regular mean curvature flow with a forcing term g 
that satisfies 

g(.,t) G W 3 '°°(R d ), d t g E W 1 ' 00 ^ x (0,T)). (13) 
Given e > 0, let u e be solution of ([7]) : 

8 t u = An - ^W'{u) + -^2W{u)g, 

and let d£l e {t) = {i£ H d ; u e (x,t) = Assume that the potential W is given 
by W(s) = ^(1 — s 2 ) 2 . Then there exist eo > and a constant C depending only 
on T such that for all e in (0, eo], the following estimate holds: 

V* G [0,T], 8n e (t) C |x G R d ; dist(x, < Ce 2 |loge| 2 } . (14) 
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Notations and assumptions. Let T > 0. For all t in [0, T], let Q(t) be a mean 
curvature flow with a forcing term g that satisfies (|13p . In the sequel, we will for 
convenience identify the signed distance to Q(t) to a function d: ~R d x [0, T] — > R 
defined by 



d(x,i) = d(x,Q(t)) 



'dist(ar,n(t)) if x e R d \ O(i) , 
ifxe<9Q(i), 
k -dist(x,fi(i)) ifxeil(t). 



We assume that dQ(t) is smooth enough so that d satisfies 

d, d t d, d t d xx d e C°(A), (15) 

where A is a tubular neighborhood of d£l(t). We assume that d£l(t) is oriented by 
the outward normal vector n defined at a point x of t?0(t) by n(x,t) = Vd(x,t). 
We denote by «i, . . . , k^-i the principal curvatures of d£l(t), and we set 

d-l d-l 

n(x,t) = Kj(x, t), h(x,t) = i<j{x, t). 
i=l i=l 

We choose k to be negative for convex balls. The evolution of d£l(t) is defined 
by V n (x,t) = n(x,t) + g(x,t) for all (x,t) in dQ(t) x [0, T], where 14 denote the 
normal velocity. 

Given D > 0, we define a tubular neighborhood A(i) of dd(t) by 

A(t) = {x G R d ; |d(x,i)| < £>} , (16) 

and we set 

A= |J A(i)x{i}. 
te[o,T\ 

If D is sufficiently small, one can associate to any point (x, t) of A a unique 
projection s(x,t) on dil(t) such that 

dist(s(x, t), x) = \d(x,t)\. 

For any scalar or vector function / defined on dCl(t), we denote by / its extension 
on A, defined by f(x,t) = f (s(x,t),t) . If / is real- valued, then we clearly have 
Vd • V/ = on A. It follows from CED that 

INIl°°(A), ||^||ioc( A ), ||V^||ioo (A ), ||Afe[| L00 ( A ) < +oo. (17) 
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Moreover, geometric properties of the distance function d imply 
d-l 



Ad(x, t) = Y j - — M Ki ^-] = -k(x, t) - d(x, t)h(x, t) + 0(d(x, t) 2 ) , 
1 — d{x, t)Ki(x, t) ' 

d t d(x, t) = -V n (x, t) = -R(x, t) - g(x, t). 



These estimates show that the motion of d£l(t) can be described by an equation 
on d inside the whole A (see [2]): 

V(x,t) G A, d t d(x,t) - Ad(x,t) = -g(x,t) + d(x,t)h(x,t) +0(d{x,t) 2 ). (18) 

We denote by q the profile function associated with the double well potential 

W: 



C(0) = 



q = argmin j^Qc' 2 + ^(C)) ds ; ( G ^ oc (R), JunJ = =Fl, 

The Euler equation for this problem writes q" = W'(q). More precisely, as W is 
smooth, q is strictly decreasing and we have q' = — y / 2W(q). If W is defined by 
W(s) = ^(1 — s 2 ) 2 , q is given by q(s) = — tanh(s). In this case, there exists a 
positive constant C such that \q — 1| < —Cq'. 
Let £ in H 2 (R) be solution of equation 

H"(s)-W"(q(s))as)=sq'(s) (19) 

with initial condition £(0) = 0. Existence and uniqueness results for this equation 
may be found in section 3 of [5], along with the following estimate: 

\e(s)\<-C(l + s 2 W(s). (20) 

Comparison lemma. Our proof of convergence relies on the following lemma: 

Lemma 1. Let e > 0, and let u and v in L 2 (0, T; H 2 (R d )) n H 1 (0, T; L 2 (R d )) 
be such that 

d t u -Au + -^W'{u) - -^2W(u)g > d t v - Av + -^W{v) - ^ y / 2W(v)g (21) 

in R d x (0,T), and u(x,0) > v(x,0) for x in K d . Then u > v in K d x (0,T). 

Proof. Let e = max(t> — u, 0). Multiplying (|2ip by e and integrating over R rf , we 
get 

jM,t)\\l^) <^(W'(u)-W'(v),e) LHnd) 



dt 

2 



y/2W(u) - y/2W(v] 



9, e 



e W / / L 2 (R d ) 

<^(W' g>e (u)-W' g>e (v),e) L2{Rdy 
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where W g e is defined by 



W'gjs) = W'{a) - e^2W(sjg, W(0) = 0. 
The idea is then to decompose W ge under the form 

W' =W' T +W' T 

where W eL is Lipschitz continuous on R, and W e 7 is nondecreasing. More 



precisely, when W(s) = ^(1 — s 2 ) 2 , we can use 

W^ L (s) = W^(s)X[-i,i](s) and W g>eJ {s) = W g>e (s)(l - X[-i,i]W). 

which satisfy the previous assumption if ||<7||l°° < §■ Then, noticing that e(x, 0) 
by assumption, we obtain 



l e (->*)lll,2( R d) < ~2 

e ./o 



'W' g ^u)-W' g ^v),e{.,r)) 



L 2 (R d ) 



(It 



< 



sup{lip« e>i )} f\\e{., 

;eR d JO 



e 2 



r )IL2(Rd) dr - 



Note that the sup is bounded just as ||<7||l°°- Gronwall's lemma implies that for 
almost every t in (0,T), ||e(., i) 1 1 ^2 (R, d ) = 0> an d e = almost everywhere in 
R d x(0,T). □ 

Construction of a subsolution. Using our previous asymptotic expansion of 
u e , we now build a subsolution to problem ([7j). Let 5 > 3 be a fixed integer. For 
all e > 0, we set s t = 5|log e|. Since <j(s) = — tanh(s), we have 



-1 + 



2e : 



2tf 



■1 + 0(e 25 ), g'( Se ) = -(1 - q(s e ) 2 ) = 0(e 25 ), 



1 + e 28 

and it follows from (|20p that 

|£( Se )| = 0(e 25 |loge| 2 ), = 0(e 25 |loge| 2 ) 

We define two auxiliary functions q e and £ e by 

if < s < s f , 



Qe(s) 



P q {s) ifs € <s<2s e 
-1 if s > 2s e , 

—q e (—s) if s < 0, 
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and 



f (s) if < s < s t , 

Pf(a) ifs e <s<2s e 

if s > 2s e , 

-&(-«) ifs<0, 



where P 9 and Pg are polynomials of degree 3 defined in such a way that q e and 
£ e are in C 1 (R). It follows that 

\\P q + l|U«»(/ e ) + a e ||^IU-(i.) + s?||i^'|U=o(/ e ) < C(|5(a e ) + 1| + s e |g'(s e )|) , 
II^IU'»(/«) + s e |l J €lU'»(/«) + s ?ll J € , IU'»(/ e ) ^ C , (|g(s e )| + s e |g / (s e )| 
with J e = [s t , 2se]. Then we easily check that 

Ik, - g||L~(R) = ^e 25 - 1 ), ||6 - eiUoo (R) = o{e 25 - 1 ), 

together with 

q , :-W'(q e ) = o(e 2S - 1 ), q' e + y/M{£) = o^" 1 ), (22) 

and 

C-W"fe)ee-^ = o(e 2<5 - 1 ). 
For e > 0, we introduce the modified distance function c£~ defined by: 

V(x,i) G R d x [0,T], d^(x,t) = d{x,t) + ci (t)e 2 |log e| 2 , 

where ci is a positive continuous function, independent of e, that will be deter- 
mined later. For t in [0, T], we introduce the sets 

A~(t) = jx G H d ; |<C"0M)l < 2-5e|loge|} 

and 

A e "= (J A e -(t)x{f}. 

ie[0,T] 

It is then possible to find eo > depending only on 5, ci, D, such that 

Ve<e , ViG[0,r], A^(i) C A(i), (23) 
where A(i) is the tubular neighborhood defined in (|16p . In particular, we see that 
V(x,t)GA7, d(x,t) = 0(e|loge|). 
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Noticing that Vd~ = Vd and Vd~ ■ Vh = in A~, it follows from JTH]) that 

- Ad~ = d t d -Ad + c[e 2 \log e\ 2 

= -9 + d~h + (4 - cJi)e 2 \\og e\ 2 + 0(e 2 |loge| 2 ). (24) 

Setting y = we define v~ on R d x [0, T] by 

g e (l/) + e 2 (h + Vd • V 5 )£e(y) - c 2 e 3 |loge| 2 in A", 
— 1 — C2e 3 |loge| 2 in {d~ > 2<5e|log e|}, 

k +l - c 2 e 3 |loge| 2 in < -25e|log e|}, 

where c 2 is a constant independent of e that we will be determined later. In view of 
(fT5|) . we easily check that v~ belongs to L 2 (0, T; i?/ oc (R d )) n i? 1 (0, T; L 2 oc (R d )) . 
Our goal is to show that v~ is a subsolution of ([7]). 
Let w e be solution of (|7|). We will first prove that 

Vx 6 R d , v~(x,0) <u t (x,0). (25) 

To this end, we introduce w e defined by 

w e = q(y) + e 2 (h + Vd ■ Vg)£(y) - |e 3 |log e| 2 , 

and we note that when e is sufficiently small, 

v e (x,0) < We (x,0) -ye 3 |loge| 2 + o(e 2<5 ~ 1 ) 
< w e (x,0), 

so that (|25p follows from showing that 

w e (x,0) - u e (x,0) = q(y(x,Q)) - q ^ X ^ 

+ e 2 (h(x, 0) + Vd(x, 0) • Vg(x, 0)) ^(x, 0)) - |e 3 |log e| 
is non-positive. We define for convenience 

h = e 2 (S(x, 0) + Vd(x, 0) • Vg{x, 0))({y(x, 0)) . 
The following lemma is proved in section 6 of [5] (recall that q' is negative). 
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Lemma 2. Let z = } and let y = dc ^'°' > = z + ci(0)e|log e| 2 . Then for e 
sufficiently small, 



2q'(y) < q'(s) < l - q \y) 



for all s in [z,y]. 

This lemma implies that 



/, _ i d 7&°) \ _ Jd-(x,0)-ci(py\\oge\ 



i Y^)V l(0)e|loge|2) 



2 

and using (fT9|) , it follows that 

J 2 < -K C e 2 (l + y(x,0) 2 )g'(y(x,0)), 

where K = \\h(., 0) || z,o° (A(o)) + l|V<7(->0)||z,°°(Rd)- We then distinguish two cases. 
If \y\ > |loge|, then (1 + y 2 )\q'(y)\ < 0(e 2 |loge| 2 ) and L 2 is controlled by the 
negative term —^e^\\oge\ 2 , so that 

w e {x,0) - u e {x,0) = I 1+ I 2 - ye 3 |loge| 2 
</2-|6 3 |log6| 2 

<0(6 4 |l0 ge | 2 )-| e 3 |l0g6| 2 . 

If \y\ < | log e | , then L 2 is controlled by I\, and 
w e (x,0) - u e (x,0) = h + J 2 - y£ 3 |loge| 2 

< q '(^^j (|ci(o) e |io ge | 2 - o(6 2 |io ge | 2 ; 

Thus, choosing ci(0) and c 2 sufficiently large, we get the desired estimate ([25]) . 
Let us now check that 

d t v~ - Av~ + \w\v~ ) - - J2W(v7)g < (26) 



e 2 • - • e 



in R d x (0,T). 
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Case 1: (x,t) £ A e . In this case, (fT7|) implies that 

dtv7 = - e q' e (y)dtd; + e 2 (d t (h + Vd ■ Vgj)Uv) + <h + Vd ■ Vg)C{y)d t d~, 
Vi>7 = \q'MVd: + e 2 (v(h + Vd ■ Vg))Uv) + < h + ■ V^^jVC 

A«r = -^'(y) + -^(y)AC + (/* + Vd • Vg)C(y) + O(e). 

Using these equalities with we get 
5t« e - - A^- = -\q'i{y) + V e (l/)(^dT " AC) ~(h + Vd- Vg)g(y) + 0(e) 

= --^q"(y) + -q' e (y) {-a + - e 2 |ioge| 2 (ci/i - c[)j 

- (h + Vd-Vg)C{y) +0(e\\oge\ 2 ) 

= ^q'I(y) - \q'e{y)9 + (q'e(y)yh -{h + Vd- Vg)C{y)) 

- e |loge| 2 (^(y)( Cl ^ - ci)) + 0(e|loge| 2 ). 

Since (x, t) is in A~, we obtain the following estimates on the terms of order in 
425): 

g(x, t) = g(s(x, t) + dVd, tj 

= g + dVd-Vg + 0{d 2 ) 

= g + eyVd- ■ Vg - Cl e 2 |log e\ 2 Vd~ • Vg + 0(e 2 |log e| 2 ), 
\w\v-) = \W(q £ ) + W"(q e )(h + Vd • Vg)& - W" (q e )c 2 e\loge\ 2 + O(e), 

-^2W(v e )g = -^2W(q e )g + 0(e) 

= -^2W(q e )(g + eyVd t ■ Vg - Cl e 2 \log e\ 2 Vd~ ■ V<?) 
+ 0(e|loge| 2 ) 

= -~ e q' t {g + eyVd~ - Vg - Cl e 2 \\oge\ 2 Vd: - Vg) + 0{e\\oge\ 2 ) 
= -- e q' t (g + eyVd- - V g - Cl e 2 \\oge\ 2 Vd ■ Vg) + 0(e|log e| 2 ). 

Here, we used and the fact that for (x,t) in A~, 

yq' e Vd-(Vg-Vg) = 0(e\loge\). 
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Summing the previous equalities, we obtain 

8 t v e - Av e + ^W'(v e ) - ±y/2W(v £ )g = h + h + h + h + 0(e|loge| 2 ), 

with 

h = -±(q':(y)-W'( qe (y)))=o(e 2S - 3 ), 
h = - e g' t (g - 9) = 0, 

I 3 = -(h + Vd- Vg) (g(y) - W"{q e % - yq> € ) = o(e 2S - 3 ), 
h = e\\oge\ 2 (-q' e ( Cl h - c[ + Cl Vd • Vg) - c 2 W"(q e )) • 

We now determine the function c\ and the constant c 2 so that I4 is sufficiently 
negative to compensate the term of order e|loge| 2 . Letting K = \\h\\ L ao^ d ^ + 

l|Vfli|i,°o(R«ix(o,:r))> we set 

ci(t) = cexp((l + 

so that 

-e|loge|V e (ci^-t4 + Vd- Vgc^j < cie|loge| 2 ^. 

We thus have 

h < -e\\oge\ 2 c 2 (--q' e + W"{q t ) 



C2 

Noticing that —csq' e + W"(q e ) is uniformly positive for C3 large enough, we can 
choose c and c 2 such that 

dtv~ - Av: + - (W'{y-) - gyj2W(vZ)) < 

in A~. 

Case 2: (x,t) £ A~ . Here, the function v~ is given by v e = ±1 — c 2 e 3 |loge| 2 , 
which implies dtv~ = 0, Av~ = 0, 

^W'(v-) = ^{W'(±l) - W'(±l)c 2 e 3 |loge| 2 + 0(e 3 )) 
= -W'(±l)c 2 e|loge| 2 + 0(e), 
-^2W(v7)g = \ ( V2W(±l)g - (vW)'(±l) 5 c 2 e 3 |log e| 2 + 0(e 3 )) 
= 0(e 2 |loge| 2 ). 
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Noticing that W"(±l) > and choosing c 2 large enough guarantees that in 
R d \A r , 

1 



d t v~ - Av~ + -(W'{v-) - g^2W{v7)) < 0. 

To conclude, we apply the comparison principle of lemma Q] and discover 

V(x,i) G R d x [0,T], v e (x, t) < u e (x, t). 

A similar argument can be applied to 

r q e (y) + e 2 (h-Vd- Vg)Uv) + c 2 e 3 |loge| 2 in A+ 
-1 + c 2 e 3 |loge| 2 in {d+ > 2(5e|loge|}, 

+1 + c 2 e 3 |loge| 2 in {d+ < -2<5e|loge|}, 

with y = — and 

d£{x,t) = d(x,t) - ci(t)e 2 |log e| 2 , 
to show that vf is a supersolution to ([7]): 

V(x,t) G H d x [0,T], u+(x,i) > u e (x,t). 
Proof of the theorem. 

Proof of theorem [3 We choose eo so that (|23[) holds. Let t in [0, T] and x in 
dQ e (t) be given. We first show that x is in A(i). Indeed, u e (x,t) = and 

«rOM) < u e (x,i) = < v+(x,t). (27) 

Assume that x g" A(i). As A^(t) C A(t), we have x G" A^(t), and thus, for e 
sufficiently small, we deduce that v~(x,t) and vf(x,i) have the same sign. This 
contradicts (p7|) . and we conclude that x G A(t). We then notice that 

v 7 (x,t) = J^±)+O(e*)<0 



because 

e 2 ( / l _ Vd . V5)e£(y) = C , (e 2 )) 

and hence 

J^M)<o( f 2 ). 



As </(0) = —1, we get that 
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which shows that d(x,t) > 0(e | log e| ). In a similar way, noticing that 



u +(x,t)=g e (^Y^)+O(O>0, 



we get that 

< 0(e2)) 

e 

and d(x,t) < 0(e 2 |log e| 2 ). We conclude that 

\d(x,t)\ <0(e 2 |loge| 2 ), 
and ([T3]) is proved. □ 

3 Application to mean curvature flow with conserva- 
tion of the volume 

In this section, we compare two phase field models for the approximation of 
motion by mean curvature with conservation of the volume: 



f Kda. (28) 
Jdn(t) 



V n = K 

Jan(t) 

As explained in the introduction, this motion is usually approximated by the 
following phase field equation (see [6]): 

d t u = Au- -=W'(u) + \-f W'(u) dx. (29) 
e e JQ 

The last term in this equation can be understood as a Lagrange multiplier for 
the mass constraint 

A r 

u dx = 0. 



-I 



(Note that in this section, the potential W we consider has its wells at and 1.) 
In the sequel, we compare this equation to 



dt u = Au-\w\u) + \ ^Hp ^-f W\u)dx. (30) 
e e 2 J Rd y / 2W(u) dx Jr* 

derived along the same lines as ([7]). The form of the last term is again related 
to conservation of mass, since the volume average of the right-hand side is easily 
seen to vanish. 
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There is no general proof of convergence of solutions of (|29p and (|30p to the 
motion (|28p . However, (|29|) is commonly used in computations. The numerical 
experiments presented further show that (|30p conserves volume with a higher 
degree of accuracy than (I29p . Our aim in this section is to try so substantiate 
this claim, although our arguments are formal. 

In both cases, the last term could be interpreted as a forcing term, by setting 



ge{t) = — iw'(u)dx 
tew Jo 



£CW JQ 

in the first model and 



e J Rd ^/2W{u)dx 

in the second. Formally, one recovers the expressions of ([2]) and ([7]). However 
the forcing terms here depend on the solutions of (|29p and (|30p . Assuming that 
one can generalize the results of section 12.11 (notwithstanding this dependence of 
g e and g e ), we expect solutions u e and u e of ([29]) and ([30]) to have the following 
asymptotic behavior: 

n e (x,t) = gl 1 + 1 + 0{e ), (31) 

n e (x,t) = g V eW M +e 2 /i(x,t)d ^ 6W M +0(e 3 ), (32) 



where O e (t) (resp. Cl e (t)) denotes the set contained inside the level line {u e (x, t) = 
^} (resp. {u e (x,t) = i}), and q, r], £ are the profiles defined in ([3]), (|3J) and 
(fT9l) . We note that these profiles only depend on the choice of the potential W. 
Following (fT2|) . we see that as g e does not depend on x, only h appears in the 
term of order 2 of u e . 

We first establish the connection between the mass Jq u t dx (respectively 
§QU e dx) and the volume |fi e (*)| (respectively |f2 e (t)|). 

Proposition 1. Let E be a regular bounded domain of~R d , and let 

'd(x,E)' 



v e (x) 

Assume that q is symmetric, i.e. q(s) = l — q(—s), and that q decays exponentially 
to as s — > +oo. Then 



\E\ = [ v e dx + 0(e 2 ). 
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Proof. Using the co-area formula. 

'd(x,E) 



[ v e dx = [ q( - 



dx 



j h(s)ds + j h(s)(q[~) -ljds + j h{s)q\~) ds 
\E\ - h(-s)q(^Jds + j h(s)q(^J ds 



\E\ + 
\E\+e 



/ {h(s)-h(-s))c 
Jo 



ds 



+00 



(h(se) — h(—se))q(s) ds 



where h(s) = \Dx{d(x,E)<s}\(R- d ) is the perimeter of the level line s of the signed 
distance function to E. Since E is smooth, one can estimate h{se) — h{—se) = 
2se/i'(0) + 0(s 2 e 2 ) for s in (0, |loge|). Furthermore, since q is exponentially de- 
creasing to as s — > +00, all the moments J s>0 s n q(s) ds are finite. Thus, we can 
estimate 



loge| 



(h(se) — h(—se))q(s) ds 



< 



f g l (2seh'{0) + Cs 2 e 2 )q(s)ds 
Jo 



O(e). 



Moreover, since h(s) ~ s d 1 as s —* +00, and since h is bounded on (— oo,0), it 
is easy to check that 

p+00 p+00 

/ h(se)q(s) ds < Ce^ 1 I s d - l q{s) ds = Ofe*" 1 ), 

»/|loge| ./|loge| 

r+00 r+00 

/ h(-se)q(s)ds <C q(s)ds = 0(e). 



It follows that 



/ 

jR d 



v t dx = \E\ + 0(e 2 ). 



□ 



The result of proposition [T] still holds on a fixed bounded set Q that strictly 
contains E when e is sufficiently small. This again is a consequence of the expo- 
nential decay of q. Recalling the asymptotic form of u e , it follows from the above 
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proposition that, for the classical model ([29]) . 

In general, the term of order e does not vanish, since 
when g is symmetric, and so 

u e dx = p e (t)\ + 0(e) 



/ 

JQ 



This explains why we cannot expect the model (|29|) to converge to the motion 
(|28p with a better rate than 0(e). As for the model (|30|) . we have 



is 

/ u e cix = |Q e (i)| +0(e 2 ), 



which presents a higher degree of accuracy on volume conservation. 

We proved in the last section that solutions of ([7]) converge as e —* to motion 
by mean curvature with a forcing term dl|). Formally, the phase field equation 
([301) can be rewritten 



8 t u = An - 1 (V'(n) - e v / 2W>)^ • 



The following property shows that, under the assumption (|32p . g e converges to 
Kticr, which is formally consistent to the limiting motion ([IJ). 

Proposition 2. Let E be a regular bounded domain ofH d , and let 

v e (x) = q\ — 

Assume that q is symmetric, i.e. q(s) = l — q(—s), and that q decays exponentially 
to as s — > +oo. Then 



f Rd s/2W(v £ )dx JdE 



1 f Rd W'(v t )dx 
ef Kd ^/2W(^)dx 



i Kda + 0(e 2 ) 

JdE 
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Proof. To prove the first equality, recall that q satisfies y / 2W(q) = —q', and that 
q' is even. Let h: R — ► R be a continuous function, differentiable at s = 0, 
which grows polynomially in s. Since J R g'(s)ds = —1, arguing as in the proof 
of proposition [TJ it follows that 



= ~ f q'(s)h(> 
Jr 



(se) ds 



r+oo 

/ (h(se) + h(-se))q'(s)ds 
Jo 

[ 6 (h(se) + h(-se))q'(s)ds + 0(e 2 ) 
Jo 

3|l0g£| (2h(0) + Cs 2 e 2 )q'{s) ds + 0(e 2 ) 



/i(0)+O(e 2 ). 



Next, the co-area formula yields 



i/ rj9 v2w -ILiL^Mi^)) ds - 

Since E is smooth, and since the forcing term g is bounded, the function 

h : s i— * / g da 

Jd(x,E)=s 

is continuous, differentiable at s = and has polynomial growth at infinity: 
h(s) ~ when s — ► +oo. We can then apply the previous estimate to obtain 

-f g v / 2W(v e )dx= [ gda + 0(e 2 )= [ gda + 0{e 2 ). 

e Jn d Jd(x,E)=0 JdE 

We notice that the same argument with g = 1 leads to 



-I ^2W{v e )dx= [ da + 0(e 2 ) = \8E\ + 0(e 

e Jn, d Jd(x,E)=0 

so that combined with the previous equality, we obtain 
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Let us now prove the second equality. Recall that 

'd(x,EY 



\Vd(x,E), 



As q" = W'(q), it follows that 



i j ^ W'{v e ) dx = f d (^W(v e ) - AvA dx 



-L 

[([ Ad(x,E)da]-q'(-)ds. 

J-R\Jd(x,E)=s J e V e / 



The function 

s h li(s) = / Ad(x,E)da 

Jd(x,E)=s 

is not continuous on R, but it is constant on a sufficiently small neighborhood 
of (depending only on the topology of E) and grows polynomially like 
Arguing as in the first part of the proof, we obtain 

4 [ W'(v £ ) dx = h(0) + 0(e 2 ) = - I k da + 0(e 2 ), 
e jRd J dE 



and 

/ Kda + 0(e 2 ), 

JdE 



1 f Rd W'(v e )dx 



which completes the proof. □ 

Remark 2. The first correcting term £ vanishes at infinity in the expansion of 
u e , and so would higher order terms. If we assume that ([32]) holds, a more careful 
analysis based on proposition [H would show that 



1 f Rd W'(u e ) dx 
e f nd ^/2W{u e )dx 



i nda + 0{e 2 ). 



This heuristically justifies the use of (|30|) as an approximation to the motion 
(EH). 
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Remark 3. We can generalize our previous argument to the case of interfaces 
moving with normal velocity 

V n = K + g- f (k + g)da. (33) 
Jon 

The usual phase field approximation of such motions is based on the equation 

d t u = An - (w'{u) - ec w g^j + j \ W'{u) - ec w g^j dx, (34) 

where the last term can be understood as a Lagrange multiplier. As explained 
above, one cannot expect that this model should converge to the motion (|33p with 
a better rate than 0(e). Generalizing our previous analysis, we may instead con- 
sider the following modified phase field model, which should improve the accuracy: 



d t u = Au-^ (w'(u) - egs/W(u)\ 



1 



e 2 f Rd dx 

4 Numerical method and simulations 



j \ W'(u) - eg v / 2W(u)^j dx. (35) 



In this section, we describe the numerical method we use for solving 

1 

u(x,0) = uq(x), x 6 Q, 



d t u = Au - ^F(u) x 6 Q C K d , t G [0, T], 



(36) 



where F takes one of the following forms: 





= W'(u) 


- tc w g, 




= W'{u) 


- eg v / 2W(u) 




= W'(u) 


- tc w g - f 
JQ 




= W'{u) 


- egy/2W(u) 






y/2W(u) 



dx, 



j \W\u) - eg^[W{i 



u) I dx. 



$ Q ^/2W(u)dx.l Q 

The first form corresponds to the Allen-Cahn equation with a forcing term g. The 
second form corresponds to the modified approximation introduced in section [2j 
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Forms 3 and 4 are the respective forms when the volume is conserved (see (|34|) 
and ([35]) ). We assume that 

uo = ql J, 

where is a smooth bounded set of R rf strictly contained in the fixed box 
Q = [-§, If, with d = 2 or 3. We assume also that during the evolution the sets 
0, e (t) remain within Q, so that we may impose periodic boundary conditions on 
dQ to the solutions of ([36]) . 



4.1 Numerical scheme 



Equation (|36|) is numerically approximated via a splitting method between the 
diffusion and reaction terms. We take advantage of the periodicity to treat the 
diffusion part of the operator in the Fourier space. More precisely, the value 
u t (x, t n ) at time t n = to + n ^t is approximated by 



[ x > tn) = ^2 u e>p(*") exp(2z7rp • x). 



max pt|<P 



In a first step, we set 

" P ^x,t ri + ^J= ^ u e , p (t n + exp(2i7rp • x) 



max |pi-|<P 



with 

Ur 



P (tn + 7j] = u £iP (t n )exp(-47r 2 At\p\ 2 ). 



We then add the reaction term: 

uf(x, t n + l)= uf (tn + - ^-F (uf (t n + i 

In practice, the first step is performed via a fast Fourier transform, with a com- 
putational cost of 0(P d log-P). The corresponding numerical scheme turns out to 
be L°°-stable for the standard Allen-Calm equation with no forcing term, under 
the condition 

St < Me 2 , 

where M = (sup tG [ ,i] W"(t)} 1 . It can be shown that this condition is also 
sufficient for the modified potential W e<g . We impose this constraint in the fol- 
lowing computations for all the choices of F. We use the double well potential 
W(s) = \s 2 {\ — s) 2 , and P represents the number of Fourier modes in each 
dimension. 
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4.2 Numerical tests 



Convergence test with no forcing term. This test illustrates the conver- 
gence of our numerical scheme when we consider the equation 

d t u = An - -sW'Cu), 

with no forcing term, nor volume conservation. The initial set Qo is taken as 
a circle of radius Rq = 0.25. It should evolve as a circle, with radius R(t) = 
o 2 — 2t, that decreases to a point at the extinction time t ext = \Rq- Figure 
[2] represents Q,(t) at different times, for the choice of parameters P = 2 8 , At = 
l/P 2 and e = 2/ P. Fi gure [3] shows the error between calculated and theoretical 
extinction times for different values of e, in logarithmic scale. The error behaves 
like 0(e 2 |loge| 2 ) as expected. This indicates that, with this choice of parameters, 
the error due to our numerical scheme is negligible compared to the 'modeling' 
error due to the approximation of the motion by the phase field equation. 

0.0305328 

0.4 
0.3 
0.2 
0.1 


-0.1 

-0.2 
-0.3 
-0.4 
-0.5 

-0.4 -0.2 0.2 0.4 

Figure 2: Mean curvature flow of a circle: level set {u e (x,t) = ^} for different 
times. 

Convergence test with a constant forcing term. Here we compare the two 
phase field models (|2|) and (|7|) as approximations to the motion ([I]). Theoretically, 
both give an approximation order of 0(e 2 |log e| 2 ). We compare the numerical 
solutions in the simple case where the forcing term is a constant: g = C g . The 
initial condition is a circle of radius Rq. During the evolution, Q(t) also 
remains circular, and its radius R satisfies 

dR_ 1 
~dt-~R + Ca - 
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Figure 3: Mean curvature flow of a circle: error on the extinction time for different 
values of e (logarithmic scale). 

Assuming that C g < I/Rq, Q(t) decreases to a point, with extinction at the time 

text = ~ {-FT Ml " C g Ro) + R 

We represent on figure H] the error on the extinction time for different values of 
e, in logarithmic scale. We choose C g = 2 and C g = —2 respectively. Both 
models give comparable results, and as expected by the theory, we again observe 
a 0(e 2 |loge| 2 ) error. 



le-05 
0.001 



/ 






0( £ ) 




0(„ 2 ) 




0(£ J |lo ge | 2 ) 




classical — i 




modified — * — 




Figure 4: Mean curvature flow of a circle with a constant forcing term C g : error 
on the extinction time for different values of e (logarithmic scale). Left: C g = 2. 
Right: C g = -2. 



Conservation of the volume with no forcing term. Here the initial con- 
figuration r?o is the union of two disjoints circles of respective radii ro and Ro, 
with ro < Rq. As it evolves by conserved mean curvature flow (|28p . 0, remains 
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the union of two circles, with radii r and R solutions of 



(dr _ 1 2 

dt r r + R ' 



1 2 

R + r + R' 



It is easy to check that the smallest circle decreases and disappears at extinction 
time 



rpRp | Rj + rj L l J 1 | 2r i? 



2 ' 4 "V" ' (^o-r ) 2 , 
Meanwhile, the radius of the initially larger circle grows to a maximal value 



+ R 



o 



at extinction time. We presents results for ro = 0.1, Rq = 0.15, t ext = 0.0133, 
and for the choice of numerical parameters P = 2 8 , At = 2~ 16 . The evolution 
of r and R is plotted on figure [5] for both models ([29]) and ([30]) and for different 
choices of e:. Figure [6] depicts the error on extinction time in logarithmic scale. 
The graph clearly shows that the error on extinction time is of order e for the 



classical model, while it scales like e for the modified model ([30 



Figure 5: Conserved mean curvature flow of two disjoint circles: evolution of the 
radii against time, for different values of e. Left: classical model (129p . Right: 
modified model ([30]) . 



Conservation of the volume with a non-zero forcing term. Volume losses 
may become important when approximating forced mean curvature motion with 
the classical phase field model (|29|) . The purpose of this test is to illustrate this 
point. We choose g to be an isotropic forcing term: g(x) = c g cos{9>ir\x\) . The 
initial configuration $7q is the circle of radius Rq = 0.25 centered at 0. It should 
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Figure 6: Conserved mean curvature flow of two disjoint circles: error on the 
extinction time against e for the two models (|29p and (|30p . 

remain stationary (i.e. Q(t) = Qq for all t) whatever the value of the constant c g . 
Figure[7]represents the computed evolutions using respectively fl34f) and (|35p . The 
numerical parameters are P = 2 8 , e = 2/P and At = 1/P 2 . Clearly, the value of 
c g has a significant impact on the results when using (|34l) . Comparatively, the 
choice of c g as a negligible impact on the evolutions computed with ([35]) . This 
confirms the arguments developed in section [3j 
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Figure 7: Conserved mean curvature flow of a circle with an additional isotropic 
forcing term g(x) = c g cos(8tt\x\): stationary shape for different choices of c g . 
Left: classical model (fM|) . Right: modified model ([35]) . 

An example in 3D. Here we illustrate the benefits of our approach on a 
classical three-dimensional example: the evolution of a torus with conservation 
of the volume and no additional forcing term. This example provides a good 
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test case: because of the high values taken by the mean curvature, standard 
approaches may fail to reproduce the motion correctly. One also need to handle 
the topological change when we move from a toric shape to a spherical one. 

We clearly observe on figure [8] that the classical model (I34p leads to significant 
volume losses compared to our modified model (I35p . We plot on figure [9] the 
volume against time for both approaches. The volume error goes up to 30% for 
the classical model, whereas it is always strictly below 5% for ours. We notice 
that, in both cases, the error decreases in the second part of the evolution. Indeed, 
it is clear that the numerical error is maximal when the average mean curvature 
is maximal; when the topological change occurs, the average mean curvature 
instantly jumps to a smaller value, as the points where the mean curvature is the 
highest just disappear from the surface. 



o 
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Figure 8: Evolution of a torus by mean curvature flow with conservation of the 
volume. First line: classical model ([29]) . Second line: modified model ([30]) . 
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5 Conclusion 



We introduced in this article a modified phase field model for the approximation 
of mean curvature flow with a forcing term. We rigorously proved its convergence 
with the same order as the classical Allen-Cahn equation: 0(e 2 |log e| 2 ). 

We formally derived this model to the case of conserved mean curvature 
flow. We observed numerically an 0(e 2 ) error for the conservation of the vol- 
ume, whereas the classical conserved Allen-Cahn equation just showed an 0(e) 
error in our simulations. 

Acknowledgements. The authors would like to thank Eric Bonnetier and 
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